###############################################################################################################################################
######################################### Concentration Complete with Daily estimated Betas ###################################################
###############################################################################################################################################

#Three digit NAICS Code:
load("HHI_3.RData")

#Four digit NAICS Code:
load("HHI_4.RData")


#### Daily
load("Daily_Beta_Down.RData")
load("Daily_Beta_M.RData")
load("Daily_Beta_Up.RData")

#With Weighting:
load("Weighting_Daily_Beta_M.RData")
load("Weighting_Daily_Beta_Up.RData")
load("Weighting_Daily_Beta_Down.RData")

#Shrinkage without Weighting:
load("Shrinking_Beta_M_Daily.RData")
load("Shrinking_Beta_Up_Daily.RData")
load("Shrinking_Beta_Down_Daily.RData")

#Shrinkage with Weighting:
load("Weighting_Shrinking_Beta_M_Daily.RData")
load("Weighting_Shrinking_Beta_Up_Daily.RData")
load("Weighting_Shrinking_Beta_Down_Daily.RData")

#With Weighting and Different Windows
load("Weighting_Daily_Beta_M_6M.RData")
load("Weighting_Daily_Beta_Up_6M.RData")
load("Weighting_Daily_Beta_Down_6M.RData")

load("Weighting_Daily_Beta_M_24M.RData")
load("Weighting_Daily_Beta_Up_24M.RData")
load("Weighting_Daily_Beta_Down_24M.RData")

load("Weighting_Daily_Beta_M_60M.RData")
load("Weighting_Daily_Beta_Up_60M.RData")
load("Weighting_Daily_Beta_Down_60M.RData")

load("I_Vol.RData")

######
# this is the Link Table from the CRSP-Compustat Merged database
# it contains the elements: gvkey,LPERMNO,LINKDT,LINKENDDT
Link_Table = read.table("LinkTable.txt", header=T, sep=",")
#####

#######
# this is a Compustat file containing the following elements:
# gvkey,datadate,fyear,indfmt,consol,popsrc,datafmt,tic,cusip,conm,curcd,sale,exchg,cik,costat,naicsh,naics,sic
Compustat = read.csv("0a4a9373fde91d3a.csv", header=T)

#The NAICS Code changes over time, so we need to match this as well with a matrix that is as large as the CRSP Beta (...)
load("PERMNO_NAICS3.RData")


library(data.table)
library(stringi)

setDT(HHI_3)
#Delete all HHI_3 with NA in the NAICS_3
HHI_3 = HHI_3[!is.na(HHI_3$NAICS_3),]
HHI_3 = HHI_3[!is.na(HHI_3$V1),]

HHI_3 = HHI_3[HHI_3$fyear>1971,]
#All NAICS Values per year, we have around 27 to 29 industries in each tercile:
#Low HHI Threshold:
HHI_3_Low  = HHI_3[,NAICS_3[V1<quantile(V1, probs = c(0.33),na.rm=T)],by="fyear"]
#HHI_3[,sum(V1<quantile(V1, probs = c(0.33),na.rm=T)),by="fyear"]
#High HHI Threshold:
HHI_3_High = HHI_3[,NAICS_3[V1>quantile(V1, probs = c(0.67),na.rm=T)],by="fyear"]
#Daily
Beta_M_Daily = apply(as.matrix(Beta_M_Daily,rownames = NULL),2,as.numeric)
Beta_Up_Daily = apply(as.matrix(Beta_Up_Daily,rownames = NULL),2,as.numeric)
Beta_Down_Daily = apply(as.matrix(Beta_Down_Daily,rownames = NULL),2,as.numeric)
#Weighting
Weighting_Beta_M_Daily = apply(as.matrix(Weighting_Beta_M_Daily,rownames = NULL),2,as.numeric)
Weighting_Beta_Up_Daily = apply(as.matrix(Weighting_Beta_Up_Daily,rownames = NULL),2,as.numeric)
Weighting_Beta_Down_Daily = apply(as.matrix(Weighting_Beta_Down_Daily,rownames = NULL),2,as.numeric)
#Shrinkage:
Shrinking_Beta_M_Daily = apply(as.matrix(setDT(Shrinking_Beta_M_Daily),rownames = NULL),2,as.numeric)
Shrinking_Beta_Up_Daily = apply(as.matrix(setDT(Shrinking_Beta_Up_Daily),rownames = NULL),2,as.numeric)
Shrinking_Beta_Down_Daily = apply(as.matrix(setDT(Shrinking_Beta_Down_Daily),rownames = NULL),2,as.numeric)
#Shrinkage with Weighting:
Weighting_Shrinking_Beta_M_Daily = apply(as.matrix(setDT(Weighting_Shrinking_Beta_M_Daily),rownames = NULL),2,as.numeric)
Weighting_Shrinking_Beta_Up_Daily = apply(as.matrix(setDT(Weighting_Shrinking_Beta_Up_Daily),rownames = NULL),2,as.numeric)
Weighting_Shrinking_Beta_Down_Daily = apply(as.matrix(setDT(Weighting_Shrinking_Beta_Down_Daily),rownames = NULL),2,as.numeric)
#Weighting with different time frames:
Weighting_Beta_M_Daily_6M = apply(as.matrix(Weighting_Beta_M_Daily_6M,rownames = NULL),2,as.numeric)
Weighting_Beta_Up_Daily_6M = apply(as.matrix(Weighting_Beta_Up_Daily_6M,rownames = NULL),2,as.numeric)
Weighting_Beta_Down_Daily_6M = apply(as.matrix(Weighting_Beta_Down_Daily_6M,rownames = NULL),2,as.numeric)
Weighting_Beta_M_Daily_24M = apply(as.matrix(Weighting_Beta_M_Daily_24M,rownames = NULL),2,as.numeric)
Weighting_Beta_Up_Daily_24M = apply(as.matrix(Weighting_Beta_Up_Daily_24M,rownames = NULL),2,as.numeric)
Weighting_Beta_Down_Daily_24M = apply(as.matrix(Weighting_Beta_Down_Daily_24M,rownames = NULL),2,as.numeric)
Weighting_Beta_M_Daily_60M = apply(as.matrix(Weighting_Beta_M_Daily_60M,rownames = NULL),2,as.numeric)
Weighting_Beta_Up_Daily_60M = apply(as.matrix(Weighting_Beta_Up_Daily_60M,rownames = NULL),2,as.numeric)
Weighting_Beta_Down_Daily_60M = apply(as.matrix(Weighting_Beta_Down_Daily_60M,rownames = NULL),2,as.numeric)



#Delete all non-matching PERMNOs:
#replace the Name format by NA
colnames(Beta_M_Daily)[1] = "NaN"
colnames(Beta_Up_Daily)[1] = "NaN"
colnames(Beta_Down_Daily)[1] = "NaN"

colnames(Weighting_Beta_M_Daily)[1] = "NaN"
colnames(Weighting_Beta_Up_Daily)[1] = "NaN"
colnames(Weighting_Beta_Down_Daily)[1] = "NaN"

colnames(Shrinking_Beta_M_Daily)[1] = "NaN"
colnames(Shrinking_Beta_Up_Daily)[1] = "NaN"
colnames(Shrinking_Beta_Down_Daily)[1] = "NaN"

colnames(Weighting_Shrinking_Beta_M_Daily)[1] = "NaN"
colnames(Weighting_Shrinking_Beta_Up_Daily)[1] = "NaN"
colnames(Weighting_Shrinking_Beta_Down_Daily)[1] = "NaN"

colnames(Weighting_Beta_M_Daily_6M)[1] = "NaN"
colnames(Weighting_Beta_Up_Daily_6M)[1] = "NaN"
colnames(Weighting_Beta_Down_Daily_6M)[1] = "NaN"

colnames(Weighting_Beta_M_Daily_24M)[1] = "NaN"
colnames(Weighting_Beta_Up_Daily_24M)[1] = "NaN"
colnames(Weighting_Beta_Down_Daily_24M)[1] = "NaN"

colnames(Weighting_Beta_M_Daily_60M)[1] = "NaN"
colnames(Weighting_Beta_Up_Daily_60M)[1] = "NaN"
colnames(Weighting_Beta_Down_Daily_60M)[1] = "NaN"
#that do not appear in PERMNO_NAICS3 but appear in Beta_M_Daily
Beta_M_Daily = Beta_M_Daily[,as.numeric(colnames(Beta_M_Daily))%in%PERMNO_NAICS3[1,]]
Beta_Up_Daily = Beta_Up_Daily[,as.numeric(colnames(Beta_Up_Daily))%in%PERMNO_NAICS3[1,]]
Beta_Down_Daily = Beta_Down_Daily[,as.numeric(colnames(Beta_Down_Daily))%in%PERMNO_NAICS3[1,]]

Weighting_Beta_M_Daily = Weighting_Beta_M_Daily[,as.numeric(colnames(Weighting_Beta_M_Daily))%in%PERMNO_NAICS3[1,]]
Weighting_Beta_Up_Daily = Weighting_Beta_Up_Daily[,as.numeric(colnames(Weighting_Beta_Up_Daily))%in%PERMNO_NAICS3[1,]]
Weighting_Beta_Down_Daily = Weighting_Beta_Down_Daily[,as.numeric(colnames(Weighting_Beta_Down_Daily))%in%PERMNO_NAICS3[1,]]

Shrinking_Beta_M_Daily = Shrinking_Beta_M_Daily[,as.numeric(colnames(Shrinking_Beta_M_Daily))%in%PERMNO_NAICS3[1,]]
Shrinking_Beta_Up_Daily = Shrinking_Beta_Up_Daily[,as.numeric(colnames(Shrinking_Beta_Up_Daily))%in%PERMNO_NAICS3[1,]]
Shrinking_Beta_Down_Daily = Shrinking_Beta_Down_Daily[,as.numeric(colnames(Shrinking_Beta_Down_Daily))%in%PERMNO_NAICS3[1,]]

Weighting_Shrinking_Beta_M_Daily = Weighting_Shrinking_Beta_M_Daily[,as.numeric(colnames(Weighting_Shrinking_Beta_M_Daily))%in%PERMNO_NAICS3[1,]]
Weighting_Shrinking_Beta_Up_Daily = Weighting_Shrinking_Beta_Up_Daily[,as.numeric(colnames(Weighting_Shrinking_Beta_Up_Daily))%in%PERMNO_NAICS3[1,]]
Weighting_Shrinking_Beta_Down_Daily = Weighting_Shrinking_Beta_Down_Daily[,as.numeric(colnames(Weighting_Shrinking_Beta_Down_Daily))%in%PERMNO_NAICS3[1,]]

Weighting_Beta_M_Daily_6M = Weighting_Beta_M_Daily_6M[,as.numeric(colnames(Weighting_Beta_M_Daily_6M))%in%PERMNO_NAICS3[1,]]
Weighting_Beta_Up_Daily_6M = Weighting_Beta_Up_Daily_6M[,as.numeric(colnames(Weighting_Beta_Up_Daily_6M))%in%PERMNO_NAICS3[1,]]
Weighting_Beta_Down_Daily_6M = Weighting_Beta_Down_Daily_6M[,as.numeric(colnames(Weighting_Beta_Down_Daily_6M))%in%PERMNO_NAICS3[1,]]

Weighting_Beta_M_Daily_24M = Weighting_Beta_M_Daily_24M[,as.numeric(colnames(Weighting_Beta_M_Daily_24M))%in%PERMNO_NAICS3[1,]]
Weighting_Beta_Up_Daily_24M = Weighting_Beta_Up_Daily_24M[,as.numeric(colnames(Weighting_Beta_Up_Daily_24M))%in%PERMNO_NAICS3[1,]]
Weighting_Beta_Down_Daily_24M = Weighting_Beta_Down_Daily_24M[,as.numeric(colnames(Weighting_Beta_Down_Daily_24M))%in%PERMNO_NAICS3[1,]]

Weighting_Beta_M_Daily_60M = Weighting_Beta_M_Daily_60M[,as.numeric(colnames(Weighting_Beta_M_Daily_60M))%in%PERMNO_NAICS3[1,]]
Weighting_Beta_Up_Daily_60M = Weighting_Beta_Up_Daily_60M[,as.numeric(colnames(Weighting_Beta_Up_Daily_60M))%in%PERMNO_NAICS3[1,]]
Weighting_Beta_Down_Daily_60M = Weighting_Beta_Down_Daily_60M[,as.numeric(colnames(Weighting_Beta_Down_Daily_60M))%in%PERMNO_NAICS3[1,]]

#that do not appear in Beta_M_Daily but in PERMNO_NAICS3
PERMNO_NAICS3 = PERMNO_NAICS3[,PERMNO_NAICS3[1,]%in%as.numeric(colnames(Beta_M_Daily))]
PERMNO = PERMNO_NAICS3[1,][-1]
#same thing now for I_Vol
I_Vol = I_Vol[,I_Vol[1,]%in%PERMNO_NAICS3[1,]]

PERMNO_NAICS3 = PERMNO_NAICS3[PERMNO_NAICS3[,1]>19720101,]
Beta_M = Beta_M_Daily[Beta_M_Daily[,1]>197112,]
Beta_Up = Beta_Up_Daily[Beta_Up_Daily[,1]>197112,]
Beta_Down = Beta_Down_Daily[Beta_Down_Daily[,1]>197112,]

Weighting_Beta_M = Weighting_Beta_M_Daily[Weighting_Beta_M_Daily[,1]>197112,]
Weighting_Beta_Up = Weighting_Beta_Up_Daily[Weighting_Beta_Up_Daily[,1]>197112,]
Weighting_Beta_Down = Weighting_Beta_Down_Daily[Weighting_Beta_Down_Daily[,1]>197112,]

Shrinking_Beta_M = Shrinking_Beta_M_Daily[Shrinking_Beta_M_Daily[,1]>197112,]
Shrinking_Beta_Up = Shrinking_Beta_Up_Daily[Shrinking_Beta_Up_Daily[,1]>197112,]
Shrinking_Beta_Down = Shrinking_Beta_Down_Daily[Shrinking_Beta_Down_Daily[,1]>197112,]

Weighting_Shrinking_Beta_M = Weighting_Shrinking_Beta_M_Daily[Weighting_Shrinking_Beta_M_Daily[,1]>197112,]
Weighting_Shrinking_Beta_Up = Weighting_Shrinking_Beta_Up_Daily[Weighting_Shrinking_Beta_Up_Daily[,1]>197112,]
Weighting_Shrinking_Beta_Down = Weighting_Shrinking_Beta_Down_Daily[Weighting_Shrinking_Beta_Down_Daily[,1]>197112,]

Weighting_Beta_M_6M = Weighting_Beta_M_Daily_6M[Weighting_Beta_M_Daily_6M[,1]>197112,]
Weighting_Beta_Up_6M = Weighting_Beta_Up_Daily_6M[Weighting_Beta_Up_Daily_6M[,1]>197112,]
Weighting_Beta_Down_6M = Weighting_Beta_Down_Daily_6M[Weighting_Beta_Down_Daily_6M[,1]>197112,]

Weighting_Beta_M_24M = Weighting_Beta_M_Daily_24M[Weighting_Beta_M_Daily_24M[,1]>197112,]
Weighting_Beta_Up_24M = Weighting_Beta_Up_Daily_24M[Weighting_Beta_Up_Daily_24M[,1]>197112,]
Weighting_Beta_Down_24M = Weighting_Beta_Down_Daily_24M[Weighting_Beta_Down_Daily_24M[,1]>197112,]

Weighting_Beta_M_60M = Weighting_Beta_M_Daily_60M[Weighting_Beta_M_Daily_60M[,1]>197112,]
Weighting_Beta_Up_60M = Weighting_Beta_Up_Daily_60M[Weighting_Beta_Up_Daily_60M[,1]>197112,]
Weighting_Beta_Down_60M = Weighting_Beta_Down_Daily_60M[Weighting_Beta_Down_Daily_60M[,1]>197112,]

I_Vol = I_Vol[I_Vol[,1]>19720101,]

#Remove the first row of NA's:
#PERMNO_NAICS3 = PERMNO_NAICS3[-1,]

library(sandwich)
library(lmtest)

#Panel Regression; Build the Panel in the following way:
#time in yyyymm   Industry Code (NAICS)   HHI(annually), fill NA's    Beta_M    Beta_CF   Beta_DR   

#Time, add this one for every PERMNO (which is every column in the matrix)
PANEL_NAICS = data.frame(Time = rep(PERMNO_NAICS3[2:dim(PERMNO_NAICS3)[1],1],dim(PERMNO_NAICS3)[2]-1), NAICS = c(PERMNO_NAICS3[2:dim(PERMNO_NAICS3)[1],-1]), PERMNO = rep(PERMNO,each=dim(PERMNO_NAICS3)[1]-1),
                         HHI = rep(NA,dim(PERMNO_NAICS3)[2]-1), Beta_M = c(Beta_M[,-1]), 
                         Beta_Down = c(Beta_Down[,-1]), Beta_Up = c(Beta_Up[,-1]), Weighting_Beta_M = c(Weighting_Beta_M[,-1]), Weighting_Beta_Up = c(Weighting_Beta_Up[,-1]), Weighting_Beta_Down = c(Weighting_Beta_Down[,-1])
                         , Shrinking_Beta_M = c(Shrinking_Beta_M[,-1]), Shrinking_Beta_Up = c(Shrinking_Beta_Up[,-1]), Shrinking_Beta_Down = c(Shrinking_Beta_Down[,-1])
                         , Weighting_Shrinking_Beta_M = c(Weighting_Shrinking_Beta_M[,-1]), Weighting_Shrinking_Beta_Up = c(Weighting_Shrinking_Beta_Up[,-1]), Weighting_Shrinking_Beta_Down = c(Weighting_Shrinking_Beta_Down[,-1])
                         , Weighting_Beta_M_6M = c(Weighting_Beta_M_6M[,-1]), Weighting_Beta_Up_6M = c(Weighting_Beta_Up_6M[,-1]), Weighting_Beta_Down_6M = c(Weighting_Beta_Down_6M[,-1])
                         , Weighting_Beta_M_24M = c(Weighting_Beta_M_24M[,-1]), Weighting_Beta_Up_24M = c(Weighting_Beta_Up_24M[,-1]), Weighting_Beta_Down_24M = c(Weighting_Beta_Down_24M[,-1]), Weighting_Beta_M_60M = c(Weighting_Beta_M_60M[,-1]), Weighting_Beta_Up_60M = c(Weighting_Beta_Up_60M[,-1]), Weighting_Beta_Down_60M = c(Weighting_Beta_Down_60M[,-1]), I_Vol = c(I_Vol[2:dim(I_Vol)[1],-1]))

PANEL_NAICS = PANEL_NAICS[!is.na(PANEL_NAICS$NAICS),]
PANEL_NAICS = PANEL_NAICS[order(PANEL_NAICS$Time),]

HHI_3 = HHI_3[HHI_3$fyear<2019,]

options(warn=2)
for(i in 1:nrow(HHI_3)){
  Temp = PANEL_NAICS[HHI_3$NAICS_3[i] == PANEL_NAICS$NAICS,]
  Temp_HHI = HHI_3[HHI_3$NAICS_3[i] == HHI_3$NAICS_3,]
  Temp_HHI = Temp_HHI[order(Temp_HHI$fyear),]
  Expand_fyear = data.frame(fyear = rep(NA,length(unique(as.numeric(unlist(lapply(strsplit(as.character(Temp$Time),split=""),function(x) stri_paste(head(x,n=4),collapse=""))))))))
  Expand_fyear$fyear = unique(as.numeric(unlist(lapply(strsplit(as.character(Temp$Time),split=""),function(x) stri_paste(head(x,n=4),collapse="")))))
  
  Temp_HHI = merge(Temp_HHI,Expand_fyear,by="fyear",all.x=T,all.y=T)
  #Fill NA's from HHI!! By the Date:
  library(stringi)
  #How often does each year occur:
  Temp_Second = cbind(Temp_HHI$V1,table(as.numeric(unlist(lapply(strsplit(as.character(Temp$Time),split=""),function(x) stri_paste(head(x,n=4),collapse=""))))))
  Temp_Second = unname(rep(Temp_Second[,1],Temp_Second[,2]))
  
  Temp$HHI = Temp_Second
  
  PANEL_NAICS[HHI_3$NAICS_3[i] == PANEL_NAICS$NAICS,] = Temp
  print(paste("iteration",i,"of",nrow(HHI_3)))
}
options(warn=1)

save(PANEL_NAICS, file = "PANEL_NAICS_Daily.RData")
